Brain_regions_downstream
Rasmus Rydbirk
2026-01-12
Setup
library(magrittr)
library(dplyr)
library(ggplot2)
library(cowplot)
library(ggrepel)
library(ggpubr)
library(scHelper)
library(biomaRt)
library(pheatmap)
library(SummarizedExperiment)
library(qs)
system("sudo apt-get update\nsudo apt-get install libnetcdf-dev -y")
library(DEP)
tt <- Sys.time()Load data
tmp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")
dat.clean.orig <- tmp$dat.clean
df.clean <- tmp$expDesign.clean
universe <- tmp$universeMetadata
# Metadata
metadata <- read.delim("/work/proteomics/02_data/Metadata_Cohort_100723.csv", sep = ";", dec = ",", header = T) %>%
mutate(intern = gsub("K", "", .$K.no.),
Diagnosis = gsub("Unknown", "MSA", .$Diagnosis)) # 1418 should be MSA
dat.clean.orig@colData$subtype <- metadata$Subtype[match(dat.clean.orig@colData$intern, metadata$intern)]REGION-SPECIFIC PROTEINS WITH IMPUTATION
Clean data
imp <- dat.clean.orig@elementMetadata@listData %>%
as.data.frame() %>%
mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43
dat.tmp <- dat.clean.orig %>%
.[rownames(.) %in% rownames(imp)]Calculate DEPs
# dat.tmp <- dat.clean
dat.tmp@elementMetadata$ID <- dat.tmp@elementMetadata$CER_ID
dat.tmp@elementMetadata$name <- dat.tmp@elementMetadata$CER_name
all_contrasts <- test_diff(dat.tmp, type = "all")## Tested contrasts: CER_vs_FC, CER_vs_MED, CER_vs_STR, CER_vs_SN, FC_vs_MED, FC_vs_STR, FC_vs_SN, MED_vs_STR, MED_vs_SN, STR_vs_SNdep_all <- add_rejections(all_contrasts, alpha = 0.05, lfc = 2)Number of significant proteins
comb <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
comb %>%
lapply(\(x) dep_all@elementMetadata %>%
as.data.frame() %>%
.[, paste(x, "significant", sep = "_")]) %>%
setNames(comb) %>%
sapply(sum) %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(line = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(x = "", y = "No. significant proteins", fill = "Comparison") +
scale_fill_manual(values = ggsci::pal_cosmic()(10)) +
geom_text(aes(y = value + 8, label = value))
PCA
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, 45)) +
ggplot2::guides(shape = "none")## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
Correlation
plot_cor(dep_all, significant = TRUE, lower = 0.6, upper = 1, pal = "YlOrRd")
Clustering of DEPs per sample
ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE,
k = 3, col_limit = 3, show_row_names = FALSE,
indicate = c("condition"))
protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
rowIds <- ComplexHeatmap::row_order(ht) %>%
lapply(\(x) protIds[x])
go.res <- rowIds %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")## ## # go.res %>%
# names() %>%
# lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
# No significant terms, instead we ouput top hits
go.res %>%
lapply(getElement, "result") %>%
lapply(dplyr::slice, seq(10))## $`3`
## ID
## GO:0090263 GO:0090263
## GO:0060828 GO:0060828
## GO:0030111 GO:0030111
## GO:0003018 GO:0003018
## GO:0071260 GO:0071260
## GO:0071496 GO:0071496
## GO:0060070 GO:0060070
## GO:0090288 GO:0090288
## GO:0051882 GO:0051882
## GO:0006972 GO:0006972
## Description
## GO:0090263 positive regulation of canonical Wnt signaling pathway
## GO:0060828 regulation of canonical Wnt signaling pathway
## GO:0030111 regulation of Wnt signaling pathway
## GO:0003018 vascular process in circulatory system
## GO:0071260 cellular response to mechanical stimulus
## GO:0071496 cellular response to external stimulus
## GO:0060070 canonical Wnt signaling pathway
## GO:0090288 negative regulation of cellular response to growth factor stimulus
## GO:0051882 mitochondrial depolarization
## GO:0006972 hyperosmotic response
## GeneRatio BgRatio RichFactor FoldEnrichment zScore pvalue
## GO:0090263 6/191 46/6858 0.13043478 4.683360 4.242328 0.001597978
## GO:0060828 9/191 100/6858 0.09000000 3.231518 3.804608 0.001785780
## GO:0030111 10/191 123/6858 0.08130081 2.919167 3.635090 0.002172139
## GO:0003018 11/191 146/6858 0.07534247 2.705228 3.524939 0.002427273
## GO:0071260 5/191 35/6858 0.14285714 5.129394 4.145255 0.002616427
## GO:0071496 5/191 35/6858 0.14285714 5.129394 4.145255 0.002616427
## GO:0060070 9/191 107/6858 0.08411215 3.020111 3.564524 0.002849182
## GO:0090288 5/191 36/6858 0.13888889 4.986911 4.059294 0.002970923
## GO:0051882 3/191 12/6858 0.25000000 8.976440 4.680581 0.003885627
## GO:0006972 4/191 25/6858 0.16000000 5.744921 4.022646 0.004661548
## p.adjust qvalue
## GO:0090263 0.8816214 0.8816214
## GO:0060828 0.8816214 0.8816214
## GO:0030111 0.8816214 0.8816214
## GO:0003018 0.8816214 0.8816214
## GO:0071260 0.8816214 0.8816214
## GO:0071496 0.8816214 0.8816214
## GO:0060070 0.8816214 0.8816214
## GO:0090288 0.8816214 0.8816214
## GO:0051882 0.9169262 0.9169262
## GO:0006972 0.9169262 0.9169262
## geneID
## GO:0090263 GPC5/PPM1B/PPM1A/USP8/USP47/TBL1XR1
## GO:0060828 GPC5/DKK3/PPM1B/PPM1A/RBX1/USP8/USP47/HECW1/TBL1XR1
## GO:0030111 GPC5/DKK3/PPM1B/PPM1A/RBX1/USP8/RACK1/USP47/HECW1/TBL1XR1
## GO:0003018 SLC4A3/SLC44A1/SOD3/PRKG1/SLC38A2/FAAH/MANF/GCLC/AKAP12/SLC2A1/TJP2
## GO:0071260 FAS/SLC38A2/TMEM87A/GCLC/SLC2A1
## GO:0071496 FAS/SLC38A2/TMEM87A/GCLC/SLC2A1
## GO:0060070 GPC5/DKK3/PPM1B/PPM1A/RBX1/USP8/USP47/HECW1/TBL1XR1
## GO:0090288 HRG/CADM4/PPM1A/EPN2/CASK
## GO:0051882 ALB/GCLC/RACK1
## GO:0006972 RCSD1/SLC25A23/ICOSLG/SLC2A1
## Count
## GO:0090263 6
## GO:0060828 9
## GO:0030111 10
## GO:0003018 11
## GO:0071260 5
## GO:0071496 5
## GO:0060070 9
## GO:0090288 5
## GO:0051882 3
## GO:0006972 4
##
## $`2`
## ID Description
## GO:0006897 GO:0006897 endocytosis
## GO:0072350 GO:0072350 tricarboxylic acid metabolic process
## GO:0006886 GO:0006886 intracellular protein transport
## GO:0043523 GO:0043523 regulation of neuron apoptotic process
## GO:0021955 GO:0021955 central nervous system neuron axonogenesis
## GO:0043547 GO:0043547 positive regulation of GTPase activity
## GO:0001845 GO:0001845 phagolysosome assembly
## GO:0048169 GO:0048169 regulation of long-term neuronal synaptic plasticity
## GO:0043114 GO:0043114 regulation of vascular permeability
## GO:0060627 GO:0060627 regulation of vesicle-mediated transport
## GeneRatio BgRatio RichFactor FoldEnrichment zScore pvalue
## GO:0006897 26/180 444/6858 0.05855856 2.231081 4.403457 8.400259e-05
## GO:0072350 4/180 12/6858 0.33333333 12.700000 6.659447 1.926603e-04
## GO:0006886 24/180 437/6858 0.05491991 2.092449 3.874543 4.224967e-04
## GO:0043523 10/180 125/6858 0.08000000 3.048000 3.793681 1.576887e-03
## GO:0021955 4/180 20/6858 0.20000000 7.620000 4.867305 1.599618e-03
## GO:0043547 7/180 70/6858 0.10000000 3.810000 3.879400 2.270127e-03
## GO:0001845 3/180 11/6858 0.27272727 10.390909 5.117218 2.511748e-03
## GO:0048169 4/180 23/6858 0.17391304 6.626087 4.436916 2.749374e-03
## GO:0043114 4/180 24/6858 0.16666667 6.350000 4.310246 3.232495e-03
## GO:0060627 19/180 366/6858 0.05191257 1.977869 3.156545 3.275084e-03
## p.adjust qvalue
## GO:0006897 0.2171467 0.2154003
## GO:0072350 0.2490135 0.2470108
## GO:0006886 0.3640513 0.3611235
## GO:0043523 0.7573412 0.7512504
## GO:0021955 0.7573412 0.7512504
## GO:0043547 0.7573412 0.7512504
## GO:0001845 0.7573412 0.7512504
## GO:0048169 0.7573412 0.7512504
## GO:0043114 0.7573412 0.7512504
## GO:0060627 0.7573412 0.7512504
## geneID
## GO:0006897 SH3GL3/AP1S1/AP2S1/ACTG1/TSC2/SGIP1/YES1/ARL8B/EPS15L1/CSNK1G1/REPS1/SH3GL2/PPP3R1/PPP3CA/LYST/CORO1A/PIK3C3/SNX12/BCR/APP/RALB/RAB1A/HRAS/VAC14/PRKCG/VPS33B
## GO:0072350 IDH3A/IDH3B/IREB2/IDH2
## GO:0006886 STXBP1/AP1S1/AP2S1/TSC2/TOMM22/TOMM70/PPP3R1/PPP3CA/RAB24/IPO13/NF1/SEC24B/COPZ1/PIK3C3/XPO7/BCR/APBA1/NDEL1/TIMM13/RAB1A/STX7/USP9X/EXOC6B/VPS33B
## GO:0043523 STXBP1/PHB1/DNAJC5/CORO1A/NF1/AARS1/SARM1/APP/HRAS/PRKCG
## GO:0021955 ARHGAP35/DCC/KIFBP/NDEL1
## GO:0043547 RGS6/RALGAPB/NF1/BCR/NDEL1/HRAS/ABR
## GO:0001845 ARL8B/CORO1A/VPS33B
## GO:0048169 FXR1/NF1/APP/HRAS
## GO:0043114 ARHGAP35/YES1/SH3GL2/BCR
## GO:0060627 STXBP1/SH3GL3/AP2S1/ACTG1/TSC2/SGIP1/CHMP3/PPP3R1/PPP3CA/DNAJC5/RAB3C/PRKCB/CORO1A/SNX12/BCR/APP/VAC14/PRKCG/TBC1D4
## Count
## GO:0006897 26
## GO:0072350 4
## GO:0006886 24
## GO:0043523 10
## GO:0021955 4
## GO:0043547 7
## GO:0001845 3
## GO:0048169 4
## GO:0043114 4
## GO:0060627 19
##
## $`1`
## ID Description
## GO:0006397 GO:0006397 mRNA processing
## GO:0016071 GO:0016071 mRNA metabolic process
## GO:0042776 GO:0042776 proton motive force-driven mitochondrial ATP synthesis
## GO:0051028 GO:0051028 mRNA transport
## GO:0015986 GO:0015986 proton motive force-driven ATP synthesis
## GO:0007158 GO:0007158 neuron cell-cell adhesion
## GO:0006403 GO:0006403 RNA localization
## GO:0008380 GO:0008380 RNA splicing
## GO:0050657 GO:0050657 nucleic acid transport
## GO:0050658 GO:0050658 RNA transport
## GeneRatio BgRatio RichFactor FoldEnrichment zScore pvalue
## GO:0006397 16/128 318/6858 0.05031447 2.695755 4.270240 0.0002583117
## GO:0016071 19/128 445/6858 0.04269663 2.287605 3.873439 0.0005451515
## GO:0042776 6/128 60/6858 0.10000000 5.357813 4.675393 0.0008247422
## GO:0051028 7/128 84/6858 0.08333333 4.464844 4.406210 0.0009203638
## GO:0015986 6/128 65/6858 0.09230769 4.945673 4.407695 0.0012612538
## GO:0007158 3/128 13/6858 0.23076923 12.364183 5.655717 0.0015840344
## GO:0006403 8/128 131/6858 0.06106870 3.271947 3.620659 0.0030021369
## GO:0008380 13/128 292/6858 0.04452055 2.385327 3.336245 0.0030064054
## GO:0050657 7/128 105/6858 0.06666667 3.571875 3.662357 0.0033499474
## GO:0050658 7/128 105/6858 0.06666667 3.571875 3.662357 0.0033499474
## p.adjust qvalue
## GO:0006397 0.3686057 0.3686057
## GO:0016071 0.3686057 0.3686057
## GO:0042776 0.3686057 0.3686057
## GO:0051028 0.3686057 0.3686057
## GO:0015986 0.4041057 0.4041057
## GO:0007158 0.4229372 0.4229372
## GO:0006403 0.4482364 0.4482364
## GO:0008380 0.4482364 0.4482364
## GO:0050657 0.4482364 0.4482364
## GO:0050658 0.4482364 0.4482364
## geneID
## GO:0006397 SNRNP40/DDX42/ADAR/CSDC2/PNN/RBM25/ARVCF/SAFB/TRA2B/CSTF3/HNRNPF/SNRNP70/DDX39B/THOC6/CD2BP2/SRRM1
## GO:0016071 SNRNP40/DDX42/CNOT10/ADAR/CSDC2/PNN/RBM25/ARVCF/SAFB/TRA2B/CSTF3/CNOT9/HNRNPF/GIGYF2/SNRNP70/DDX39B/THOC6/CD2BP2/SRRM1
## GO:0042776 SDHD/NDUFB1/NDUFA11/NDUFB9/NDUFA8/ATP5MF
## GO:0051028 NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2
## GO:0015986 SDHD/NDUFB1/NDUFA11/NDUFB9/NDUFA8/ATP5MF
## GO:0007158 CNTN4/NLGN3/ASTN2
## GO:0006403 NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2/STAU2
## GO:0008380 SNRNP40/DDX42/PNN/RBM25/ARVCF/TRA2B/HNRNPF/SNRNP70/DDX39B/RTCB/THOC6/CD2BP2/SRRM1
## GO:0050657 NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2
## GO:0050658 NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2
## Count
## GO:0006397 16
## GO:0016071 19
## GO:0042776 6
## GO:0051028 7
## GO:0015986 6
## GO:0007158 3
## GO:0006403 8
## GO:0008380 13
## GO:0050657 7
## GO:0050658 7Clustering of DEPs per comparison
plot_heatmap(dep_all, type = "contrast", kmeans = TRUE,
k = 3, col_limit = 5, show_row_names = FALSE)
Volcano plots
combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)## $V1## Warning: ggrepel: 16 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V2## Warning: ggrepel: 58 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V3## Warning: ggrepel: 161 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V4## Warning: ggrepel: 137 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V5
##
## $V6
##
## $V7
##
## $V8## Warning: ggrepel: 48 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V9## Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V10
Extract results
res <- dep_all@elementMetadata@listData %>%
as.data.frame() %>%
mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x)))
# res2 <- get_results(dep_all)
res.split <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>% filter(!!sym(paste0(cc, "_p.adj")) <= 0.05, abs(!!sym(paste0(cc, "_diff"))) >= 1.5)), .)GO on DEPs
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: KPNA2,APLF,HELQ,KAT5,SMARCAD1,MEAF6## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: SLC25A36,SLC5A1,SPIDR,STAT6,IL4,RRP7BP## --> return NULL...res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: NEIL3,H1-8,SLC25A33,KDM1A,EIF6,CLCF1## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: APTX,H1-7,APLF,TERF2IP,SESN2,NDFIP1## --> return NULL...res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Get markers
all.markers <- res.split %>%
names() %>%
lapply(\(comp) {
up <- strsplit(comp, "_vs_") %>% sget(1)
down <- strsplit(comp, "_vs_") %>% sget(2)
diff <- paste0(comp, "_diff")
up.df <- res.split[[comp]] %>%
filter(!!sym(diff) > 0) %>%
mutate(area = up) %>%
dplyr::select(name, !!sym(diff), area)
down.df <- res.split[[comp]] %>%
filter(!!sym(diff) < 0) %>%
mutate(area = down) %>%
dplyr::select(name, !!sym(diff), area)
out <- rbind(up.df, down.df) %>%
dplyr::rename(l2fc = !!sym(diff)) %>%
mutate(l2fc = abs(l2fc))
return(out)
}) %>%
bind_rows()
markers.sum <- all.markers %>%
dplyr::group_by(area, name) %>%
summarize(n = n(),
mean_l2fc = mean(l2fc),
max_l2fc = max(l2fc)) %>%
data.frame()## `summarise()` has grouped output by 'area'. You can override using the
## `.groups` argument.markers.split <- markers.sum %>%
arrange(-n, -mean_l2fc) %$%
split(., area)Plot top markers
markers.split %>%
lapply(dplyr::slice, seq(10)) %>%
lapply(pull, name) %>%
Map(\(genes, name) plot_single(dep = dep_all, proteins = genes, type = "centered") + guides(col = "none") + labs(title = name), genes = ., name = names(.))## $CER
##
## $FC
##
## $MED
##
## $SN
##
## $STR
Onthology on top markers
markers.enrich <- markers.split %>%
lapply(dplyr::slice, seq(50)) %>%
lapply(pull, name) %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
markers.enrich %>%
lapply(head, 20) %>%
lapply(knitr::kable)## $CER
##
##
## | |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:----------|:----------|:------------------------------------------------------------------------------------|:---------|:--------|----------:|--------------:|---------:|--------:|---------:|---------:|:--------------------------------------------------------------------------------------------------------------------------------|-----:|
## |GO:0006397 |GO:0006397 |mRNA processing |19/49 |318/6858 | 0.0597484| 8.362341| 11.404180| 0.00e+00| 0.0000000| 0.0000000|WDR33/ALKBH5/SRSF2/SNU13/SNRPD3/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1 | 19|
## |GO:0016071 |GO:0016071 |mRNA metabolic process |21/49 |445/6858 | 0.0471910| 6.604815| 10.371323| 0.00e+00| 0.0000000| 0.0000000|WDR33/ALKBH5/SRSF2/SNU13/FUS/FTO/SNRPD3/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1 | 21|
## |GO:0000377 |GO:0000377 |RNA splicing, via transesterification reactions with bulged adenosine as nucleophile |16/49 |215/6858 | 0.0744186| 10.415567| 11.898906| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1 | 16|
## |GO:0000398 |GO:0000398 |mRNA splicing, via spliceosome |16/49 |215/6858 | 0.0744186| 10.415567| 11.898906| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1 | 16|
## |GO:0000375 |GO:0000375 |RNA splicing, via transesterification reactions |16/49 |219/6858 | 0.0730594| 10.225328| 11.769988| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1 | 16|
## |GO:1903241 |GO:1903241 |U2-type prespliceosome assembly |8/49 |22/6858 | 0.3636364| 50.894249| 19.883107| 0.00e+00| 0.0000000| 0.0000000|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF | 8|
## |GO:0008380 |GO:0008380 |RNA splicing |17/49 |292/6858 | 0.0582192| 8.148309| 10.589312| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/FUS/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1 | 17|
## |GO:0000245 |GO:0000245 |spliceosomal complex assembly |9/49 |60/6858 | 0.1500000| 20.993878| 13.194882| 0.00e+00| 0.0000000| 0.0000000|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39 | 9|
## |GO:0036260 |GO:0036260 |RNA capping |5/49 |12/6858 | 0.4166667| 58.316327| 16.856755| 0.00e+00| 0.0000010| 0.0000009|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF | 5|
## |GO:0022618 |GO:0022618 |protein-RNA complex assembly |9/49 |134/6858 | 0.0671642| 9.400244| 8.330170| 3.00e-07| 0.0000273| 0.0000249|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39 | 9|
## |GO:0000966 |GO:0000966 |RNA 5'-end processing |5/49 |23/6858 | 0.2173913| 30.425909| 11.990796| 5.00e-07| 0.0000341| 0.0000311|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF | 5|
## |GO:0071826 |GO:0071826 |protein-RNA complex organization |9/49 |141/6858 | 0.0638298| 8.933565| 8.074464| 5.00e-07| 0.0000351| 0.0000321|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39 | 9|
## |GO:0000387 |GO:0000387 |spliceosomal snRNP assembly |5/49 |28/6858 | 0.1785714| 24.992711| 10.791242| 1.30e-06| 0.0000820| 0.0000748|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF | 5|
## |GO:0001510 |GO:0001510 |RNA methylation |5/49 |31/6858 | 0.1612903| 22.574062| 10.212246| 2.20e-06| 0.0001295| 0.0001182|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF | 5|
## |GO:0045934 |GO:0045934 |negative regulation of nucleobase-containing compound metabolic process |14/49 |474/6858 | 0.0295359| 4.133816| 5.998458| 3.30e-06| 0.0001808| 0.0001650|H1-10/H1-4/HMGB1/ZNF512B/HDGF/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/SRSF7/CCAR1/KHDRBS1 | 14|
## |GO:0022613 |GO:0022613 |ribonucleoprotein complex biogenesis |10/49 |238/6858 | 0.0420168| 5.880638| 6.500691| 5.20e-06| 0.0002528| 0.0002307|SNU13/SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39 | 10|
## |GO:0051253 |GO:0051253 |negative regulation of RNA metabolic process |13/49 |423/6858 | 0.0307329| 4.301346| 5.945801| 5.30e-06| 0.0002528| 0.0002307|H1-4/HMGB1/ZNF512B/HDGF/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/SRSF7/CCAR1/KHDRBS1 | 13|
## |GO:0043484 |GO:0043484 |regulation of RNA splicing |7/49 |117/6858 | 0.0598291| 8.373626| 6.823942| 1.67e-05| 0.0007526| 0.0006869|SRSF2/FUS/SF3B3/SRSF7/SRSF3/SNRNP70/KHDRBS1 | 7|
## |GO:0043414 |GO:0043414 |macromolecule methylation |5/49 |49/6858 | 0.1020408| 14.281549| 7.914593| 2.28e-05| 0.0009724| 0.0008874|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF | 5|
## |GO:1903311 |GO:1903311 |regulation of mRNA metabolic process |8/49 |178/6858 | 0.0449438| 6.290300| 6.066323| 3.20e-05| 0.0012974| 0.0011841|ALKBH5/SRSF2/FUS/FTO/SRSF7/SRSF3/SNRNP70/KHDRBS1 | 8|
##
## $FC
##
##
## | |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:----------|:----------|:-----------------------------------------------|:---------|:--------|----------:|--------------:|--------:|---------:|---------:|---------:|:---------------------------------------------------------------------------------|-----:|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission |13/47 |495/6858 | 0.0262626| 3.832108| 5.433655| 0.0000178| 0.0076320| 0.0069256|EPHA4/CAMKV/SLC17A7/CHRM1/CAMK2A/NPY/GABRA1/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A | 13|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling |13/47 |495/6858 | 0.0262626| 3.832108| 5.433655| 0.0000178| 0.0076320| 0.0069256|EPHA4/CAMKV/SLC17A7/CHRM1/CAMK2A/NPY/GABRA1/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A | 13|
## |GO:0021954 |GO:0021954 |central nervous system neuron development |5/47 |51/6858 | 0.0980392| 14.305382| 7.922187| 0.0000226| 0.0076320| 0.0069256|PLXNA1/EPHA4/GABRB1/NPY/SLC4A10 | 5|
## |GO:0050808 |GO:0050808 |synapse organization |11/47 |394/6858 | 0.0279188| 4.073766| 5.220100| 0.0000521| 0.0132177| 0.0119943|SYNPO/PLXNA1/EPHA4/ACTN1/ICAM5/CAMKV/LINGO2/GABRA1/AMIGO1/TNC/HOMER1 | 11|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission |10/47 |346/6858 | 0.0289017| 4.217193| 5.101156| 0.0000911| 0.0157871| 0.0143259|EPHA4/CAMKV/CHRM1/CAMK2A/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A | 10|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling |10/47 |347/6858 | 0.0288184| 4.205040| 5.089615| 0.0000933| 0.0157871| 0.0143259|EPHA4/CAMKV/CHRM1/CAMK2A/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A | 10|
## |GO:0060078 |GO:0060078 |regulation of postsynaptic membrane potential |5/47 |80/6858 | 0.0625000| 9.119681| 6.067974| 0.0001990| 0.0283848| 0.0257575|GABRB1/CHRM1/GABRA1/RGS7BP/STX1A | 5|
## |GO:0060560 |GO:0060560 |developmental growth involved in morphogenesis |6/47 |129/6858 | 0.0465116| 6.786739| 5.511427| 0.0002249| 0.0283848| 0.0257575|PLXNA1/ISLR2/CPNE5/RASAL1/OLFM1/TNC | 6|
## |GO:0007186 |GO:0007186 |G protein-coupled receptor signaling pathway |9/47 |318/6858 | 0.0283019| 4.129667| 4.747154| 0.0002531| 0.0283848| 0.0257575|NECAB2/CHRM1/CAMK2A/NPY/RGS7BP/CHGA/PRKAR2B/GNG2/HOMER1 | 9|
## |GO:0021953 |GO:0021953 |central nervous system neuron differentiation |5/47 |86/6858 | 0.0581395| 8.483424| 5.800984| 0.0002797| 0.0283848| 0.0257575|PLXNA1/EPHA4/GABRB1/NPY/SLC4A10 | 5|
## |GO:0006821 |GO:0006821 |chloride transport |4/47 |56/6858 | 0.0714286| 10.422492| 5.881012| 0.0005480| 0.0459986| 0.0417410|GABRB1/SLC17A7/GABRA1/SLC4A10 | 4|
## |GO:0043270 |GO:0043270 |positive regulation of monoatomic ion transport |5/47 |101/6858 | 0.0495050| 7.223510| 5.233951| 0.0005891| 0.0459986| 0.0417410|SCN3B/CHRM1/CAMK2A/AMIGO1/HOMER1 | 5|
## |GO:1990138 |GO:1990138 |neuron projection extension |5/47 |101/6858 | 0.0495050| 7.223510| 5.233951| 0.0005891| 0.0459986| 0.0417410|PLXNA1/ISLR2/CPNE5/RASAL1/OLFM1 | 5|
##
## $MED
##
##
## |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:--|:-----------|:---------|:-------|----------:|--------------:|------:|------:|--------:|------:|:------|-----:|
##
## $SN
##
##
## | |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:----------|:----------|:-------------------------------------------|:---------|:--------|----------:|--------------:|--------:|---------:|---------:|---------:|:--------------------------------------------------------------|-----:|
## |GO:0042220 |GO:0042220 |response to cocaine |5/46 |38/6858 | 0.1315789| 19.616705| 9.456079| 0.0000046| 0.0045342| 0.0040917|SLC6A4/PPP1R1B/GAD2/HTR1B/CNR1 | 5|
## |GO:0043279 |GO:0043279 |response to alkaloid |5/46 |65/6858 | 0.0769231| 11.468227| 6.967993| 0.0000665| 0.0183148| 0.0165274|SLC6A4/PPP1R1B/GAD2/HTR1B/CNR1 | 5|
## |GO:0006837 |GO:0006837 |serotonin transport |3/46 |13/6858 | 0.2307692| 34.404682| 9.906060| 0.0000771| 0.0183148| 0.0165274|SLC6A4/HTR1B/CNR1 | 3|
## |GO:0007631 |GO:0007631 |feeding behavior |4/46 |37/6858 | 0.1081081| 16.117509| 7.576450| 0.0000995| 0.0183148| 0.0165274|PMCH/HTR1B/CNR1/CARTPT | 4|
## |GO:0003013 |GO:0003013 |circulatory system process |9/46 |292/6858 | 0.0308219| 4.595146| 5.159008| 0.0001115| 0.0183148| 0.0165274|SLC2A13/SLC6A4/SCN3B/HTR1B/ATP8A1/CNR1/KCNMA1/CARTPT/SCN2B | 9|
## |GO:0043269 |GO:0043269 |regulation of monoatomic ion transport |8/46 |227/6858 | 0.0352423| 5.254166| 5.356073| 0.0001120| 0.0183148| 0.0165274|SLC6A4/TESC/DPP6/SCN3B/HTR1B/CNR1/HPCA/SCN2B | 8|
## |GO:0015844 |GO:0015844 |monoamine transport |4/46 |45/6858 | 0.0888889| 13.252174| 6.775774| 0.0002159| 0.0289610| 0.0261347|SLC6A4/HTR1B/CNR1/CARTPT | 4|
## |GO:0007617 |GO:0007617 |mating behavior |3/46 |19/6858 | 0.1578947| 23.540046| 8.084329| 0.0002539| 0.0289610| 0.0261347|SLC6A4/PPP1R1B/CNR1 | 3|
## |GO:0023061 |GO:0023061 |signal release |8/46 |260/6858 | 0.0307692| 4.587291| 4.845692| 0.0002857| 0.0289610| 0.0261347|SLC6A4/VGF/WLS/SYT5/SLC32A1/HTR1B/CNR1/CARTPT | 8|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission |11/46 |495/6858 | 0.0222222| 3.313044| 4.389996| 0.0003247| 0.0289610| 0.0261347|PMCH/PDYN/SLC6A4/VGF/SYT5/SLC32A1/GAD2/HTR1B/CNR1/CARTPT/SCN2B | 11|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling |11/46 |495/6858 | 0.0222222| 3.313044| 4.389996| 0.0003247| 0.0289610| 0.0261347|PMCH/PDYN/SLC6A4/VGF/SYT5/SLC32A1/GAD2/HTR1B/CNR1/CARTPT/SCN2B | 11|
## |GO:0010765 |GO:0010765 |positive regulation of sodium ion transport |3/46 |24/6858 | 0.1250000| 18.635870| 7.111704| 0.0005180| 0.0410515| 0.0370452|TESC/SCN3B/SCN2B | 3|
## |GO:0019098 |GO:0019098 |reproductive behavior |3/46 |25/6858 | 0.1200000| 17.890435| 6.952064| 0.0005859| 0.0410515| 0.0370452|SLC6A4/PPP1R1B/CNR1 | 3|
## |GO:0050433 |GO:0050433 |regulation of catecholamine secretion |3/46 |25/6858 | 0.1200000| 17.890435| 6.952064| 0.0005859| 0.0410515| 0.0370452|HTR1B/CNR1/CARTPT | 3|
## |GO:0071300 |GO:0071300 |cellular response to retinoic acid |3/46 |27/6858 | 0.1111111| 16.565217| 6.658917| 0.0007381| 0.0482710| 0.0435602|TNC/SLC6A4/TESC | 3|
##
## $STR
##
##
## | |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:----------|:----------|:-------------------------------------------------------|:---------|:--------|----------:|--------------:|---------:|---------:|---------:|---------:|:------------------------------------------------------------------------------------------------------------------|-----:|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission |19/47 |495/6858 | 0.0383838| 5.600774| 8.826999| 0.0000000| 0.0000001| 0.0000001|CHRM1/CAMKV/AKAP5/SLC17A7/GABRG2/KCNIP2/ACHE/NPY/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/SLC1A2/LAMP5 | 19|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling |19/47 |495/6858 | 0.0383838| 5.600774| 8.826999| 0.0000000| 0.0000001| 0.0000001|CHRM1/CAMKV/AKAP5/SLC17A7/GABRG2/KCNIP2/ACHE/NPY/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/SLC1A2/LAMP5 | 19|
## |GO:0043270 |GO:0043270 |positive regulation of monoatomic ion transport |9/47 |101/6858 | 0.0891089| 13.002317| 10.093909| 0.0000000| 0.0000074| 0.0000060|CHRM1/AKAP5/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/THY1/SLC9A1 | 9|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission |14/47 |346/6858 | 0.0404624| 5.904071| 7.775855| 0.0000000| 0.0000094| 0.0000076|CHRM1/CAMKV/AKAP5/ACHE/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/LAMP5 | 14|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling |14/47 |347/6858 | 0.0403458| 5.887056| 7.760662| 0.0000000| 0.0000094| 0.0000076|CHRM1/CAMKV/AKAP5/ACHE/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/LAMP5 | 14|
## |GO:0051050 |GO:0051050 |positive regulation of transport |15/47 |461/6858 | 0.0325380| 4.747773| 6.920638| 0.0000002| 0.0000397| 0.0000321|CHRM1/AKAP5/RAB27B/SCN3B/KCNIP2/ACHE/ACTN2/ATP2B1/SYT1/STX1A/HOMER1/PLCB1/THY1/SLC9A1/SLC1A2 | 15|
## |GO:0043269 |GO:0043269 |regulation of monoatomic ion transport |10/47 |227/6858 | 0.0440529| 6.427969| 6.908302| 0.0000023| 0.0003719| 0.0003007|CHRM1/AKAP5/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/DPP10/THY1/SLC9A1 | 10|
## |GO:0010959 |GO:0010959 |regulation of metal ion transport |9/47 |188/6858 | 0.0478723| 6.985288| 6.912125| 0.0000040| 0.0005689| 0.0004599|AKAP5/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/DPP10/THY1/SLC9A1 | 9|
## |GO:0095500 |GO:0095500 |acetylcholine receptor signaling pathway |4/47 |22/6858 | 0.1818182| 26.529981| 9.962570| 0.0000129| 0.0016365| 0.0013230|CHRM1/ACHE/PLCB1/LY6H | 4|
## |GO:0051928 |GO:0051928 |positive regulation of calcium ion transport |5/47 |49/6858 | 0.1020408| 14.889275| 8.104878| 0.0000185| 0.0019256| 0.0015567|AKAP5/ATP2B1/HOMER1/THY1/SLC9A1 | 5|
## |GO:1905145 |GO:1905145 |cellular response to acetylcholine |4/47 |24/6858 | 0.1666667| 24.319149| 9.505860| 0.0000186| 0.0019256| 0.0015567|CHRM1/ACHE/PLCB1/LY6H | 4|
## |GO:0098926 |GO:0098926 |postsynaptic signal transduction |4/47 |25/6858 | 0.1600000| 23.346383| 9.297841| 0.0000220| 0.0019300| 0.0015603|CHRM1/ACHE/PLCB1/LY6H | 4|
## |GO:1905144 |GO:1905144 |response to acetylcholine |4/47 |25/6858 | 0.1600000| 23.346383| 9.297841| 0.0000220| 0.0019300| 0.0015603|CHRM1/ACHE/PLCB1/LY6H | 4|
## |GO:0098660 |GO:0098660 |inorganic ion transmembrane transport |12/47 |458/6858 | 0.0262009| 3.823098| 5.194924| 0.0000414| 0.0033664| 0.0027216|GABRB1/AKAP5/SLC17A7/SCN3B/GABRG2/KCNIP2/ACTN2/ATP2B1/DPP10/PLCB1/THY1/SLC9A1 | 12|
## |GO:0034764 |GO:0034764 |positive regulation of transmembrane transport |6/47 |99/6858 | 0.0606061| 8.843327| 6.529610| 0.0000517| 0.0037049| 0.0029951|AKAP5/KCNIP2/ACTN2/THY1/SLC9A1/SLC1A2 | 6|
## |GO:0050808 |GO:0050808 |synapse organization |11/47 |394/6858 | 0.0279188| 4.073766| 5.220100| 0.0000521| 0.0037049| 0.0029951|ICAM5/CAMKV/ACTN1/SYNPO/GABRG2/ACHE/ACTN2/HOMER1/PLXNA1/RAP2A/TNC | 11|
## |GO:1904062 |GO:1904062 |regulation of monoatomic cation transmembrane transport |7/47 |153/6858 | 0.0457516| 6.675845| 5.897764| 0.0000713| 0.0046683| 0.0037740|AKAP5/KCNIP2/ACTN2/HOMER1/DPP10/THY1/SLC9A1 | 7|
## |GO:1904064 |GO:1904064 |positive regulation of cation transmembrane transport |5/47 |65/6858 | 0.0769231| 11.224223| 6.879656| 0.0000738| 0.0046683| 0.0037740|AKAP5/KCNIP2/ACTN2/THY1/SLC9A1 | 5|
## |GO:0030001 |GO:0030001 |metal ion transport |11/47 |415/6858 | 0.0265060| 3.867624| 5.006252| 0.0000836| 0.0050070| 0.0040478|AKAP5/SLC17A7/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/DPP10/PLCB1/THY1/SLC9A1 | 11|
## |GO:1901699 |GO:1901699 |cellular response to nitrogen compound |10/47 |351/6858 | 0.0284900| 4.157120| 5.043880| 0.0001027| 0.0058442| 0.0047246|CHRM1/GABRB1/GABRG2/ACHE/ACTN2/ATP2B1/PLCB1/SLC9A1/SLC1A2/LY6H | 10|no.sig <- markers.enrich %>%
lapply(getElement, "result") %>%
sapply(\(x) sum(x$p.adjust <= 0.05))
markers.enrich %>%
.[no.sig > 0] %>%
names() %>%
lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
DISEASE-SPECIFIC PROTEINS WITH IMPUTATION
Prepare
dat.clean <- dat.clean.orig
dat.clean@colData@listData$area <- dat.clean@colData@listData$condition
dat.clean@colData@listData$areaRep <- dat.clean@colData@listData$replicate
dat.clean@colData@listData$condition <- dat.clean@colData@listData$Diagnosis
dat.clean@colData@listData$replicate <- dat.clean@colData@listData$diagRep
dat.clean@elementMetadata$ID <- dat.clean@elementMetadata$CER_ID
dat.clean@elementMetadata$name <- dat.clean@elementMetadata$CER_name
dat.clean@colData@rownames <- dat.clean$ID
dat.all <- dat.cleanFC
Extract results
min.l2fc = 0.1
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "FC"]
all_contrasts <- test_diff(dat.clean, type = "all")## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRLcomb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame() %>%
mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43
res.split <- comb %>%
setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
# filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc) %>%
arrange(!!sym(paste0(cc, "_p.adj")))), .)Number of significant proteins
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
Volcano plots
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)## $V1
##
## $V2
##
## $V3
##
## $V4
##
## $V5
##
## $V6
PCA
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
Clustering of DEPs per sample
plot_heatmap(dep_all, type = "centered", kmeans = TRUE,
k = 2, col_limit = 3, show_row_names = FALSE, # k 1 or 2?
indicate = c("condition"))
# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
#
# rowIds <- ComplexHeatmap::row_order(ht) %>%
# lapply(\(x) protIds[x])
#
# go.res <- rowIds %>%
# lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
#
# go.res %>%
# names() %>%
# lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))Clustering of DEPs per comparison
plot_heatmap(dep_all, type = "contrast", kmeans = TRUE,
k = 3, col_limit = 5, show_row_names = FALSE)
Onthology on top significant proteins
markers.enrich <- res.split %>%
lapply(dplyr::slice, seq(50)) %>%
.[sapply(., nrow) > 0] %>%
lapply(pull, name) %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
markers.enrich %<>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
markers.enrich %>%
names() %>%
lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))## [[1]]
##
## [[2]]
CER
Extract results
min.l2fc = 0.1
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "CER"]
all_contrasts <- test_diff(dat.clean, type = "all")## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRLcomb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame() %>%
mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43
res.split <- comb %>%
setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
# filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc) %>%
arrange(!!sym(paste0(cc, "_p.adj")))), .)Number of significant proteins
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
Volcano plots
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)## $V1
##
## $V2
##
## $V3
##
## $V4## Warning: ggrepel: 28 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V5
##
## $V6
PCA
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
Clustering of DEPs per sample
ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE,
k = 2, col_limit = 3, show_row_names = FALSE,
indicate = c("condition"))
protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
rowIds <- ComplexHeatmap::row_order(ht) %>%
lapply(\(x) protIds[x])
go.res <- rowIds %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
go.res %>%
names() %>%
lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))## [[1]]
##
## [[2]]
Clustering of DEPs per comparison
plot_heatmap(dep_all, type = "contrast", kmeans = TRUE,
k = 3, col_limit = 5, show_row_names = FALSE)
Onthology on top significant proteins
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: MIURF,XNDC1N,TP53,MRTO4,AURKC,C1QBP## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: RPS28,RPA2,TWNK,MORF4L2,OOEP,KPNA2## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: MEF2A,KPNA1,NBN,RPF2,XRCC1,IL27RA## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: RPSA,MAD2L2,MCAT,HELB,TERF2,MORF4L2## --> return NULL...res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: TBX21,VPS72,H1-8,BOP1,LIG4,RAD51## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: KMT5C,IL4,KPNA2,APLF,HELQ,KAT5## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: SMARCAD1,MEAF6,SLC25A36,SLC5A1,SPIDR,STAT6## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: IL4,SPIDR,NEIL3,H1-8,SLC25A33,KDM1A## --> return NULL...res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))## [[1]]
##
## [[2]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))## [[1]]
##
## [[2]]
STR
Extract results
min.l2fc = 2
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "STR"]
all_contrasts <- test_diff(dat.clean, type = "all")## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRLcomb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame() %>%
mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43
res.split <- comb %>%
setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj")))), .)Number of significant proteins
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
Volcano plots
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)## $V1## Warning: ggrepel: 62 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V2
##
## $V3
##
## $V4## Warning: ggrepel: 37 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V5## Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V6
PCA
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4, indicate = "condition")$data %>%
ggplot(aes(PC1, PC2, shape = condition, col = subtype)) +
geom_point(size = 3) +
theme_bw()
Clustering of DEPs per sample
ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE,
k = 2, col_limit = 3, show_row_names = F,
indicate = c("condition"))
protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
rowIds <- ComplexHeatmap::row_order(ht) %>%
lapply(\(x) protIds[x])
go.res <- rowIds %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
go.res %>%
names() %>%
lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))## [[1]]
##
## [[2]]
Clustering of DEPs per comparison
plot_heatmap(dep_all, type = "contrast", kmeans = TRUE,
k = 3, col_limit = 5, show_row_names = FALSE)
Onthology on top significant proteins
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: HDGFL2,MIURF,XNDC1N,TP53,MRTO4,AURKC## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: C1QBP,RPS28,RPA2,TWNK,MORF4L2,OOEP## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: KPNA2,MEF2A,KPNA1,NBN,RPF2,XRCC1## --> return NULL...res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: IL27RA,RPSA,MAD2L2,MCAT,HELB,TERF2## --> return NULL...res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))## [[1]]
##
## [[2]]
##
## [[3]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))## [[1]]
##
## [[2]]
STRINGdb network
Networks are shown for all evidence of interaction, not just physical evidence (stronger).
sdb <- STRINGdb::STRINGdb(species = 9606, score_threshold = 400, version = "12.0")sdb.prots <- res.split %>%
names() %>%
lapply(\(x) {
message(x)
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_"))) %>%
pull(name)
}) %>%
setNames(names(res.split)) %>%
.[sapply(., length) > 1]## PSP_vs_MSA## PSP_vs_PD## PSP_vs_CTRL## MSA_vs_PD## MSA_vs_CTRL## PD_vs_CTRLsdb.prots %>%
lapply(sdb$plot_network)
## $PSP_vs_MSA
## NULL
##
## $MSA_vs_PD
## NULL
##
## $MSA_vs_CTRL
## NULLSN
Extract results
min.l2fc = 0.5
min.padj = 0.05
dat.clean <- dat.all[, dat.all$area == "SN"]
all_contrasts <- test_diff(dat.clean, type = "all")## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRLcomb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame() %>%
mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43
res.split <- comb %>%
setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj")))), .)Number of significant proteins
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
Volcano plots
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)## $V1
##
## $V2
##
## $V3## Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##
## $V4
##
## $V5
##
## $V6
PCA
plot_pca(dep_all, x = 1, y = 2, n = 100, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
Clustering of DEPs per sample
ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE,
k = 2, col_limit = 3, show_row_names = FALSE,
indicate = c("condition"))
protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
rowIds <- ComplexHeatmap::row_order(ht) %>%
lapply(\(x) protIds[x])
go.res <- rowIds %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
go.res %>%
names() %>%
lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))## [[1]]
##
## [[2]]
Clustering of DEPs per comparison
plot_heatmap(dep_all, type = "contrast", kmeans = TRUE,
k = 2, col_limit = 5, show_row_names = FALSE)
Onthology on top significant proteins
res.go.up <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) > 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: CGAS,KMT5B,HDGFL2,MIURF,XNDC1N,TP53## --> return NULL...res.go.down <- res.split %>%
names() %>%
lapply(\(x) {
res.split[[x]] %>%
filter(!!sym(paste(x, "significant", sep = "_")),
!!sym(paste(x, "diff", sep = "_")) < 0) %>%
pull(name) %>%
clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
}) %>%
setNames(names(res.split)) %>%
.[!sapply(., is.null)] %>%
.[sapply(., nrow) > 0] %>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]## --> No gene can be mapped....## --> Expected input gene ID: MRTO4,AURKC,C1QBP,RPS28,RPA2,TWNK## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: MORF4L2,OOEP,KPNA2,MEF2A,KPNA1,NBN## --> return NULL...## --> No gene can be mapped....## --> Expected input gene ID: RPF2,XRCC1,IL27RA,RPSA,MAD2L2,MCAT## --> return NULL...res.go.up %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
res.go.down %>%
names() %>%
lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))## [[1]]
##
## [[2]]
MED
Extract results
min.l2fc = 0.1
min.padj = 0.2
dat.clean <- dat.all[, dat.all$area == "MED"]
all_contrasts <- test_diff(dat.clean, type = "all")## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRLcomb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))
dep_all <- all_contrasts
for (cc in comb) {
dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
(dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
(abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}
dep_all@elementMetadata$significant <- dep_all@elementMetadata %>%
.[, grepl("_significant", names(.))] %>%
getElement("listData") %>%
as.data.frame() %>%
apply(1, \(x) any(x))
res <- dep_all@elementMetadata@listData %>%
as.data.frame() %>%
mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43
res.split <- comb %>%
setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
arrange(!!sym(paste0(cc, "_p.adj")))), .)Number of significant proteins
res.split %>%
names() %>%
lapply(\(cc) res.split[[cc]] %>%
filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>%
setNames(names(res.split)) %>%
lapply(nrow) %>%
unlist() %>%
{data.frame(comp = names(.), value = unname(.))} %>%
ggplot(aes(comp, value, fill = comp)) +
geom_bar(stat = "identity") +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
labs(x = "", y = "No. significant proteins") +
scale_fill_manual(values = ggsci::pal_cosmic()(7))
Volcano plots
combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>%
as.data.frame() %>%
apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>%
lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)## $V1
##
## $V2
##
## $V3
##
## $V4
##
## $V5
##
## $V6
PCA
plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
scale_shape_manual(values = rep(20, ncol(dat.clean))) +
ggplot2::guides(shape = "none")## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.
Onthology on top significant proteins
markers.enrich <- res.split %>%
lapply(dplyr::slice, seq(50)) %>%
.[sapply(., nrow) > 0] %>%
lapply(pull, name) %>%
lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T)
markers.enrich %<>%
.[(sapply(., \(x) {
x@result %>%
filter(p.adjust <= 0.05) %>%
nrow()})) > 0]
markers.enrich %>%
names() %>%
lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
Session
Sys.time() - tt## Time difference of 11.93094 minssessionInfo()## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
##
## Matrix products: default
## BLAS/LAPACK: /opt/intel/oneapi/mkl/2025.3/lib/libmkl_gf_lp64.so.2; LAPACK version 3.12.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Berlin
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DEP_1.26.0 qs_0.27.3
## [3] SummarizedExperiment_1.38.1 Biobase_2.68.0
## [5] GenomicRanges_1.60.0 GenomeInfoDb_1.44.0
## [7] IRanges_2.42.0 S4Vectors_0.46.0
## [9] BiocGenerics_0.54.0 generics_0.1.4
## [11] MatrixGenerics_1.20.0 matrixStats_1.5.0
## [13] pheatmap_1.0.12 biomaRt_2.64.0
## [15] scHelper_0.0.5 ggpubr_0.6.0
## [17] ggrepel_0.9.6 cowplot_1.1.3
## [19] ggplot2_3.5.2 dplyr_1.1.4
## [21] magrittr_2.0.3
##
## loaded via a namespace (and not attached):
## [1] STRINGdb_2.20.0 bitops_1.0-9
## [3] fs_1.6.6 ProtGenerics_1.40.0
## [5] enrichplot_1.28.2 httr_1.4.7
## [7] RColorBrewer_1.1-3 doParallel_1.0.17
## [9] ggsci_3.2.0 tools_4.5.2
## [11] MSnbase_2.30.1 backports_1.5.0
## [13] R6_2.6.1 DT_0.33
## [15] lazyeval_0.2.2 GetoptLong_1.0.5
## [17] withr_3.0.2 prettyunits_1.2.0
## [19] gridExtra_2.3 preprocessCore_1.70.0
## [21] fdrtool_1.2.18 cli_3.6.5
## [23] Cairo_1.6-2 sandwich_3.1-1
## [25] labeling_0.4.3 conos_1.5.2
## [27] sass_0.4.10 mvtnorm_1.3-3
## [29] readr_2.1.5 yulab.utils_0.2.0
## [31] gson_0.1.0 DOSE_4.2.0
## [33] R.utils_2.13.0 dichromat_2.0-0.1
## [35] MetaboCoreUtils_1.16.1 plotrix_3.8-4
## [37] limma_3.64.0 rstudioapi_0.17.1
## [39] impute_1.82.0 RSQLite_2.3.11
## [41] gridGraphics_0.5-1 shape_1.4.6.1
## [43] RApiSerialize_0.1.4 gtools_3.9.5
## [45] car_3.1-3 GO.db_3.21.0
## [47] Matrix_1.7-3 MALDIquant_1.22.3
## [49] imputeLCMD_2.1 abind_1.4-8
## [51] R.methodsS3_1.8.2 lifecycle_1.0.4
## [53] yaml_2.3.10 carData_3.0-5
## [55] gplots_3.2.0 qvalue_2.40.0
## [57] SparseArray_1.8.0 BiocFileCache_2.16.0
## [59] Rtsne_0.17 grid_4.5.2
## [61] blob_1.2.4 promises_1.3.2
## [63] crayon_1.5.3 PSMatch_1.12.0
## [65] shinydashboard_0.7.3 ggtangle_0.0.6
## [67] lattice_0.22-7 mzR_2.38.0
## [69] KEGGREST_1.48.0 magick_2.8.6
## [71] pillar_1.10.2 knitr_1.50
## [73] ComplexHeatmap_2.24.0 fgsea_1.34.0
## [75] rjson_0.2.23 codetools_0.2-20
## [77] fastmatch_1.1-6 glue_1.8.0
## [79] ggfun_0.1.8 pcaMethods_2.0.0
## [81] data.table_1.17.2 MultiAssayExperiment_1.34.0
## [83] vctrs_0.6.5 png_0.1-8
## [85] treeio_1.32.0 gtable_0.3.6
## [87] gsubfn_0.7 assertthat_0.2.1
## [89] cachem_1.1.0 xfun_0.52
## [91] S4Arrays_1.8.0 mime_0.13
## [93] ncdf4_1.24 iterators_1.0.14
## [95] gmm_1.8 statmod_1.5.0
## [97] nlme_3.1-168 ggtree_3.16.0
## [99] bit64_4.6.0-1 progress_1.2.3
## [101] filelock_1.0.3 bslib_0.9.0
## [103] affyio_1.78.0 tmvtnorm_1.6
## [105] KernSmooth_2.23-26 colorspace_2.1-1
## [107] DBI_1.2.3 tidyselect_1.2.1
## [109] chron_2.3-62 bit_4.6.0
## [111] compiler_4.5.2 curl_7.0.0
## [113] httr2_1.2.1 xml2_1.4.0
## [115] DelayedArray_0.34.1 stringfish_0.16.0
## [117] caTools_1.18.3 scales_1.4.0
## [119] affy_1.86.0 rappdirs_0.3.3
## [121] stringr_1.5.1 digest_0.6.37
## [123] rmarkdown_2.29 XVector_0.48.0
## [125] htmltools_0.5.8.1 pkgconfig_2.0.3
## [127] dbplyr_2.5.1 fastmap_1.2.0
## [129] rlang_1.1.6 GlobalOptions_0.1.2
## [131] htmlwidgets_1.6.4 UCSC.utils_1.4.0
## [133] shiny_1.10.0 farver_2.1.2
## [135] jquerylib_0.1.4 zoo_1.8-14
## [137] jsonlite_2.0.0 BiocParallel_1.42.0
## [139] mzID_1.46.0 GOSemSim_2.34.0
## [141] R.oo_1.27.1 Formula_1.2-5
## [143] GenomeInfoDbData_1.2.14 ggplotify_0.1.2
## [145] patchwork_1.3.0 Rcpp_1.0.14
## [147] proto_1.0.0 ape_5.8-1
## [149] MsCoreUtils_1.20.0 sqldf_0.4-11
## [151] vsn_3.76.0 leidenAlg_1.1.5
## [153] stringi_1.8.7 MASS_7.3-65
## [155] org.Hs.eg.db_3.21.0 plyr_1.8.9
## [157] parallel_4.5.2 Biostrings_2.76.0
## [159] sccore_1.0.5 splines_4.5.2
## [161] hash_2.2.6.3 hms_1.1.3
## [163] Spectra_1.18.0 circlize_0.4.16
## [165] igraph_2.1.4 QFeatures_1.18.0
## [167] ggsignif_0.6.4 reshape2_1.4.4
## [169] XML_3.99-0.18 evaluate_1.0.3
## [171] RcppParallel_5.1.10 BiocManager_1.30.25
## [173] tzdb_0.5.0 foreach_1.5.2
## [175] httpuv_1.6.16 tidyr_1.3.1
## [177] purrr_1.0.4 clue_0.3-66
## [179] norm_1.0-11.1 BiocBaseUtils_1.10.0
## [181] broom_1.0.8 xtable_1.8-4
## [183] AnnotationFilter_1.32.0 tidytree_0.4.6
## [185] rstatix_0.7.2 later_1.4.2
## [187] tibble_3.2.1 clusterProfiler_4.16.0
## [189] aplot_0.2.5 memoise_2.0.1
## [191] AnnotationDbi_1.70.0 cluster_2.1.8.1